#ifndef __OPTIMIZER_H__
#define __OPTIMIZER_H__

#include "serialize.h"

namespace Optimizer {
	class Attribute {

	};
	
	class NurbsPatchAttribute : public Attribute {
		public:
			NurbsPatchAttribute(SNurbs* e, const std::string& strAttribute) : e(e), strAttribute(strAttribute) {}
			SNurbs* e;
			std::string strAttribute;
	};
	
	class TriangleAttribute : public Attribute {
		public:
			TriangleAttribute(STriangle* t, const std::string& strAttribute) : t(t), strAttribute(strAttribute) {}
			STriangle* t;
			std::string strAttribute;
	};

	class EdgeAttribute : public Attribute {
		public:
			EdgeAttribute(const boost::array<nfloat, 4>& x, const boost::array<nfloat, 4>& y, const boost::array<nfloat, 4>& z) : x(x), y(y), z(z) {}
			boost::array<nfloat, 4> x, y, z;
	};
	
	class ElementAttribute : public Attribute {
		public:
			ElementAttribute(Element* e, const std::string& strAttribute) : e(e), strAttribute(strAttribute) {}
			Element* e;
			std::string strAttribute;
	};

	class PrimitiveEvaluator {
		public:
			virtual nfloat f(const Attribute& attribute) const = 0; 
			virtual Vector Gradient(const Attribute& attribute) const = 0; 
			virtual Matrix Hessian(const Attribute& attribute) const = 0; 
		private:
	};

	class Evaluator {
		public:
			Evaluator(Mesh* mesh) : mesh(mesh) {}
			Evaluator(Mesh* mesh, PrimitiveEvaluator* eval) : mesh(mesh), eval(eval) {}
			virtual nfloat f(const std::string& attribute) const;
			virtual nfloat f(Node* v, const std::string& attribute) const {std::cerr << "warning: no implemented." << std::endl; return 0.0;}
			virtual Vector Gradient(const std::string& attribute) const;
			virtual SparseMatrix Hessian(const std::string& attribute) const {std::cerr << "warning: no implemented." << std::endl; SparseMatrix(0, 0, 0);}
			virtual void Initialize(const std::string& attribute) const {};
			Mesh* GetMesh() {return mesh;}
		protected:
			Mesh* mesh;
			PrimitiveEvaluator* eval;
	};
	
	class TriangleEvaluator : public Evaluator {
		public:
			TriangleEvaluator(STriangleMesh* tm, PrimitiveEvaluator* eval) : Evaluator(tm, eval) {}
			virtual nfloat f(const std::string& attribute) const;
			virtual nfloat f(Node* v, const std::string& attribute) const;
			virtual Vector Gradient(const std::string& attribute) const;
			virtual SparseMatrix Hessian(const std::string& attribute) const;
			virtual void Initialize(const std::string& attribute) const {};
		protected:
	};
	
	class EdgeEvaluator : public Evaluator {
		public:
			EdgeEvaluator(STriangleMesh* tm, PrimitiveEvaluator* eval) : Evaluator(tm, eval) {}
			virtual nfloat f(const std::string& attribute) const;
			virtual nfloat f(Node* v, const std::string& attribute) const;
			virtual Vector Gradient(const std::string& attribute) const;
			virtual SparseMatrix Hessian(const std::string& attribute) const;
			virtual void Initialize(const std::string& attribute) const {};
		protected:
			std::pair< std::vector<size_t>, boost::tuples::tuple<boost::array<nfloat, 4>, boost::array<nfloat, 4>, boost::array<nfloat, 4> > > Wedge(const SEdge* e, const std::string& attribute) const;
	};

	class NurbsSurfaceEvaluator : public Evaluator {
		public:
			NurbsSurfaceEvaluator(NurbsSurfaceMesh* nsm, PrimitiveEvaluator* eval) : Evaluator(nsm, eval) {}
			virtual nfloat f(const std::string& attribute) const;
			virtual nfloat f(Node* v, const std::string& attribute) const {std::cerr << "per vertex energy not implemented." << std::endl; return 0.0;}
			virtual Vector Gradient(const std::string& attribute) const;
			virtual SparseMatrix Hessian(const std::string& attribute) const;
			virtual void Initialize(const std::string& attribute) const {};
		protected:
	};
	
	class GradientPreconditioner {
		public:
			virtual Vector operator()(const Vector& g) const {return g;}
			virtual void Update() {}
		private:
	};

	class SobolevGradientPreconditioner : public GradientPreconditioner {
		public:
			SobolevGradientPreconditioner(Mesh* mesh, const std::string& attribute, const size_t order, const size_t maxStep = 10, const bool forward = true) : mesh(mesh), attribute(attribute), order(order) {
				const size_t n = mesh->SizeNodes();
				for(size_t i = 0; i < n; i++) {
					Node* v = static_cast<Node*>(mesh->GetNode(i));
					if(v->Get<size_t>("mark") != 0) {
						boundary.push_back(v->Id());
					}
				}
			}
			virtual Vector operator()(const Vector& g) const;
			virtual void Update();
		private:
			Mesh* mesh;
			std::string attribute;
			std::vector<size_t> boundary;
            size_t order;
			
			SparseMatrix K;
			SparseMatrix M;
			SparseMatrix Minv;
	};
	
	class SobolevGradientVariableFilterPreconditioner : public GradientPreconditioner {
		public:
			SobolevGradientVariableFilterPreconditioner(Mesh* mesh, const std::string& attribute, const size_t order, const size_t maxStep = 10, const bool forward = true) : mesh(mesh), attribute(attribute), order(order), maxStep(maxStep), forward(forward) {
				const size_t n = mesh->SizeNodes();
				for(size_t i = 0; i < n; i++) {
					Node* v = static_cast<Node*>(mesh->GetNode(i));
					if(v->Get<size_t>("mark") != 0) {
						boundary.push_back(v->Id());
					}
				}
				currentStep = 0;
			}
			virtual Vector operator()(const Vector& g) const;
		private:
			Mesh* mesh;
			std::string attribute;
			std::vector<size_t> boundary;
            size_t order;
			mutable size_t currentStep;
			size_t maxStep;
			bool forward;
	};
	
	/////////////////////Global Evaluators/////////////////////////////////////////////
	class FunctionEvaluator {
		public:
			virtual nfloat operator()(const nfloat& t) const = 0;
	};
	
	class EnergyEvaluator : public FunctionEvaluator {
		public:
			EnergyEvaluator() : mesh(0), eval(0) {}
			EnergyEvaluator(Mesh* mesh, const Vector& g, Evaluator* eval, const std::string& attribute1, const std::string& attribute2) : mesh(mesh), g(g), eval(eval), attribute1(attribute1), attribute2(attribute2) {}
			virtual nfloat operator()(const nfloat& t) const;
			void SetGradient(const Vector& g) {this->g = g;}
		private:
			Vector g;
			Evaluator* eval;
			Mesh* mesh;
			std::string attribute1, attribute2;
	};
	
	class StepSizeEstimator {
		public:
			StepSizeEstimator() {}
			StepSizeEstimator(Mesh* mesh, Evaluator* eval, const std::string& attribute1, const std::string& attribute2) {
				eeval = EnergyEvaluator(mesh, Vector(0), eval, attribute1, attribute2);
			}
			virtual std::pair<nfloat, nfloat> operator()(const nfloat t, const Vector& g);
		protected:
			EnergyEvaluator eeval;
	};
	
	class LineSearchStepSizeEstimator : public StepSizeEstimator {
		public:
			LineSearchStepSizeEstimator(Mesh* mesh, Evaluator* eval, const std::string& attribute1, const std::string& attribute2) {
				eeval = EnergyEvaluator(mesh, Vector(0), eval, attribute1, attribute2);
			}
			virtual std::pair<nfloat, nfloat> operator()(const nfloat t, const Vector& g);
	};
	
	/////////////////////Local Evaluators/////////////////////////////////////////////
	class EnergyLocalEvaluator : public FunctionEvaluator {
		public:
			EnergyLocalEvaluator() : tm(0), eval(0) {}
			EnergyLocalEvaluator(STriangleMesh* tm, const Vec3n& g, Evaluator* eval, const std::string& attribute1, const std::string& attribute2) : tm(tm), g(g), eval(eval), attribute1(attribute1), attribute2(attribute2) {}
			virtual nfloat operator()(const nfloat& t) const;
			void SetGradient(const Vec3n& g) {this->g = g;}
			void SetVertex(Node* v) {this->v = v;}
		private:
			Node* v;
			Vec3n g;
			Evaluator* eval;
			STriangleMesh* tm;
			std::string attribute1, attribute2;
	};
	
	
	class StepSizeLocalEstimator {
		public:
			StepSizeLocalEstimator() {}
			StepSizeLocalEstimator(STriangleMesh* tm, Evaluator* eval, const std::string& attribute1, const std::string& attribute2) {
				eeval = EnergyLocalEvaluator(tm, math::V3n(0.0, 0.0, 0.0), eval, attribute1, attribute2);
			}
			virtual std::pair<nfloat, nfloat> operator()(const nfloat t, const Vec3n& g, Node* v);
		protected:
			EnergyLocalEvaluator eeval;
	};
	
	class LineSearchStepSizeLocalEstimator : public StepSizeLocalEstimator {
		public:
			LineSearchStepSizeLocalEstimator(STriangleMesh* tm, Evaluator* eval, const std::string& attribute1, const std::string& attribute2) {
				eeval = EnergyLocalEvaluator(tm, math::V3n(0.0, 0.0, 0.0), eval, attribute1, attribute2);
			}
			virtual std::pair<nfloat, nfloat> operator()(const nfloat t, const Vec3n& g, Node* v);
	};

	class Config {
		public:
			Config() {}
			nfloat eps;
			nfloat epsChangeOrder;
			size_t maxIterations;
			size_t maxLoops;
			size_t minOrder;
			size_t maxOrder;
			bool flagPrecondition;
			bool flagLineSearch;
			std::string info;

			template<class Archive>
				void serialize(Archive & ar, const unsigned int version) {
					ar & boost::serialization::make_nvp("Eps", eps);
					ar & boost::serialization::make_nvp("EpsChangeOrder", epsChangeOrder);
					ar & boost::serialization::make_nvp("MaxIterations", maxIterations);
					ar & boost::serialization::make_nvp("MaxLoops", maxLoops);
					ar & boost::serialization::make_nvp("MinOrder", minOrder);
					ar & boost::serialization::make_nvp("MaxOrder", maxOrder);
					ar & boost::serialization::make_nvp("PreconditionGradient", flagPrecondition);
					ar & boost::serialization::make_nvp("LineSearch", flagLineSearch);
					ar & boost::serialization::make_nvp("Info", info);
				}

			void Init(const std::string& filename) {
				std::ifstream istr(filename.c_str());
				boost::archive::xml_iarchive ia(istr);
				using namespace boost::serialization;
				ia >> make_nvp("Config", *this);
			}

			void Save(const std::string& filename) const {
				std::ofstream ostr(filename.c_str());
				boost::archive::xml_oarchive oa(ostr);
				using namespace boost::serialization;
				oa << make_nvp("Config", *this);
			}
	};
	
	std::vector< std::pair<nfloat, nfloat> > MinimizeCG(Evaluator* eval, GradientPreconditioner* gpc, StepSizeEstimator* sse, const std::string& attribute, const nfloat eps, const nfloat epsChangeOrder, const size_t maxIterations);
	void MinimizeMultiScale(Mesh* mesh, Evaluator* eval, const Config& config);
}

#endif //__OPTIMIZER_H__
